function [ param_best, cost_best, t_eq_best ] = annealing_optimization( DATA, INFO, param_0, param_id, study_case, T_0, dT, delta_t, free, t_f1, t_f2, nit_max, nit_T, LQL, activate_vd, use_Markov)

%INFO: Simulated annealing to fit model to experimental data (Dewan et al. & Deng et al.)

%% OUT:

    % param_best: optimal parameters
    % cost_best: optimal cost
    % t_eq_best: equilibrium time with param_best

%% IN:

    % DATA: experimental data (cell)
    %    DATA(:,1): days
    %    DATA(:,2): volume
    %    DATA(:,3): standard deviation
    % INFO: treatment info cell, each row represents a curve
    %    Columns: D; t_rad; anti-CTLA-4 IT (yes:1, no:0); t_treat_c4;
    %    anti-PD-(L)1 IT (yes:1, no:0); t_treat_p1
    % param_0: (vec) initial parameters
    % param_id: (vec) optimization parameters index
    % study_case: (vec) curves index
    % T_0: (real) initial temperature
    % dT: (real) temperature delta
    % free: (vec) controls free growth
    %    free(1): (int) 1 for free behavior until an specific time or volume, 0 otherwise 
    %    free(2): (int) 1 for time, 0 for volume
    %    free(3): (real) initial volume or time
    % t_f1: (real) time allowed to reach baseline tumor volume
    % t_f2: (real) simulation time
    % nit_max: (int) number of steps in which temperature decreases
    % nit_T: (int) number of evaluations for each temperature
    % LQL: =1 if LQL, otherwise (modified) LQ (with PARAM(33)) 
    % activate_vd: =1 if vascular death (D > 15) 
    % use_Markov: =1 if Markov, otherwise deterministic

%% Algorithm

% Parameters and cost initialization
    param_op = param_0;
    param_best = param_0;
    [cost_op, t_eq_op] = least_squares(DATA, INFO, param_0, study_case, delta_t, free, t_f1, t_f2, LQL, activate_vd, use_Markov);
    cost_best = cost_op;
    t_eq_best = t_eq_op;

    T = T_0; % Initial temperature

    n = nit_max; % Number of steps in which temperature decreases
    m = nit_T;   % Number of evaluations for each temperature

    for i = 1:n         
        for j = 1:m
            % Calculate new parameters
            param_new = neighbor(param_op, param_id, T, T_0);
            % Calculate the cost associated to the new parameters
            [cost_new, t_eq] = least_squares(DATA, INFO, param_new, study_case, delta_t, free, t_f1, t_f2, LQL, activate_vd, use_Markov);
            % Calculate the survival for the new parameters
            surv = survival(cost_new, cost_op, T);
            if surv == 1
                if cost_new < cost_best
                    cost_best = cost_new;
                    param_best = param_new;
                    t_eq_best = t_eq;
                end
           
                param_op = param_new;
                cost_op = cost_new;
            end
        end %for j
    
        % cooling
        T = dT * T;

    end %for i

end


%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                               %
%     Function to obtain neighbor parameters    %
%                                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function param = neighbor(param, param_id, T, T_0)

    var = (2 * rand - 1.0) * 0.5 * sqrt(T / T_0);
    m = length(param_id);    
    id = ceil(rand * m);
    
    if param_id(id) == 34
        param(param_id(id)) = param(param_id(id)) + var * 0.05;
    else
        param(param_id(id)) = max(0, param(param_id(id)) * (1 + var));
    end

    % Values constraints

    param(3) = min(0.35, max(param(3), 0.02));              % alpha_C
    param(4) = min(param(3)/2, max(param(3)/20, param(4))); % beta_C
    param(5) = min(0.7, max(param(5), 0.03));               % phi
    param(13) = min(0.7, max(param(13), 0.03));             % sigma
    param(14) = min(5, param(14));                          % tau_1
    param(16) = min(0.5, max(1e-4, param(16)));             % alpha_T
    param(17) = param(16)/10;                               % beta_T
    param(18) = min(5, param(18) );                         % tau_2
    param(19) = min(0.7, max(param(19), 0.03));             % eta
    param(22) = min(1, max(param(22), 0.05));               % h
    param(23) = min(1e-7, param(23));                       % iota
    param(26) = max(2, min(20, param(26)));                 % r
    param(28) = min(3e-7, param(28));                       % a
    param(30) = min(1, param(30));                          % q
    %param(34) constraint: use only with the modified LQ model
    param(34) = min(0.2236, param(34));                     % beta_2
    
end


%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                               %
%     Function to asign a suvival probability to a solution     %
%                                                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function surv = survival(cost_new, cost, T)

    surv = 0;

    if cost_new < cost
        surv = 1;
    elseif rand < exp(- (cost_new - cost) / T)
        surv = 1;
    end
    
end


%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                        %
%     Objective function computation     %
%                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [cost_tot, t_eq] = least_squares(DATA, INFO, param_new, study_case, delta_t, free, t_f1, t_f2, LQL, activate_vd, use_Markov)
    
    % Solution computation   
    cost = zeros(1, length(study_case));
    par_c4 = param_new(25);
    par_p1 = param_new(35);

    for i = 1:length(study_case)
        j = study_case(i);
        D = cell2mat(INFO(j, 1));
        t_rad = cell2mat(INFO(j, 2));
        param_new(25) = cell2mat(INFO(j, 3)) * par_c4;
        param_new(35) = cell2mat(INFO(j, 5)) * par_p1;
        t_treat_c4 = cell2mat(INFO(j, 4));
        t_treat_p1 = cell2mat(INFO(j, 6));

        [sol, t_eq] = radioimmuno_response_model(param_new, delta_t, free, t_f1, t_f2, D, t_rad, t_treat_c4, t_treat_p1, LQL, activate_vd, use_Markov);

        if t_eq == -1 %Initial tumor volume not achieved
            cost_tot = 1e10;
            return
        end

        ind = round( cell2mat(DATA(j, 1)) ./ delta_t + t_eq / delta_t) + 1;
        data_y = cell2mat(DATA(j, 2));
        err = cell2mat(DATA(j, 3));

        sol = sol(ind);
        if ~isrow(data_y), data_y = data_y'; end
        if ~isrow(err), err = err'; end
        cost(i) = sum(((data_y - sol) .^ 2) ./ (err .^ 2));
    end
    
    cost_tot = sum(cost);
    
end
